	global directory ""

	cap program drop myprobit
	program myprobit
	
		args lnfj xb1 xb2
				
		quietly replace `lnfj'   = ln(normal(-`xb1')) 								if $ML_y1 == 0
		quietly replace `lnfj'   = ln(binormal(`xb1',`xb2',$rho)) 					if $ML_y1 == 2
		quietly replace `lnfj'   = ln(1-normal(-`xb1')-binormal(`xb1',`xb2',$rho))	if $ML_y1 == 1
		
	end

	cap program drop model
	program model

		global rho = `1'
	
		ml model lf myprobit (Y = residscore DoB DoB2) (residscore DoB DoB2) if DoB < 0 [aw=w10], vce(robust)
		ml maximize		

		gen xb1 = [eq1]_b[_cons] + (residscore * [eq1]_b[residscore]) 
		gen xb2 = [eq2]_b[_cons] + (residscore * [eq2]_b[residscore])	
		
		gen stay = 1 - normal(-xb1)   				if DoB <  0
		gen qual = .
			replace qual = binormal(xb1,xb2,$rho)   if DoB <  0
			replace qual = 1-normal(-xb2)			if DoB >= 0
		
		reg stay residscore							if DoB <  0, r
		global stay_model_cons  = _b[_cons] 
		global stay_model_slope = _b[residscore] 

		reg qual residscore_after residscore after, r
		global qual_model_cons   		= _b[_cons] 
		global qual_model_slope  		= _b[residscore] 
		global qual_model_after  		= _b[_cons] + _b[after] 
		global qual_model_deltaslope 	= _b[residscore] + _b[residscore_after] 	
		
		cap drop xb1 xb2 qual stay
		
	end
	
	
	*** ESTIMATE MODEL
	
	use "$directory/ROSLA_Sample" if parentsscore < . & score < . & edu16 < . & CSE_Olevel_only < ., clear
	
	global qual CSE_Olevel_only
	
	gen Y = .
		replace Y = 0	if edu16 == 0 
		replace Y = 1	if edu16 == 1 & $qual == 0 
		replace Y = 2	if edu16 == 1 & $qual == 1 
		
	cap program drop boot_model
	program define boot_model, rclass
	
		cap drop residscore
		cap drop residscore_after
	
		reg score parentsscore
		predict residscore, residuals
		
		gen residscore_after = residscore * after
	
		reg edu16 	residscore 	DoB DoB2		if DoB <  0 [aw=w10]
		local stay_actual_cons  = _b[_cons] 
		local stay_actual_slope = _b[residscore] 
		
		reg $qual residscore_after residscore after DoB DoB2 DoBafter DoB2after [aw=w10]
		local qual_actual_cons   		= _b[_cons] 
		local qual_actual_slope  		= _b[residscore] 
		local qual_actual_after  		= _b[_cons] + _b[after] 
		local qual_actual_deltaslope 	= _b[residscore] + _b[residscore_after] 					
	
		return scalar stay_actual_cons  	 = `stay_actual_cons'
		return scalar stay_actual_slope 	 = `stay_actual_slope'
		return scalar qual_actual_cons  	 = `qual_actual_cons'
		return scalar qual_actual_slope 	 = `qual_actual_slope'
		return scalar qual_actual_after 	 = `qual_actual_after'
		return scalar qual_actual_deltaslope = `qual_actual_deltaslope'
	
		*******************
		***				***
		*** RHO = -0.5  ***
		***				***
		*******************
	
		model -0.5				
		
		return scalar stay_model1_cons  	 	 = $stay_model_cons
		return scalar stay_model1_slope 	 	 = $stay_model_slope
		return scalar qual_model1_cons  	 	 = $qual_model_cons
		return scalar qual_model1_slope 	 	 = $qual_model_slope
		return scalar qual_model1_after 	 	 = $qual_model_after
		return scalar qual_model1_deltaslope  	 = $qual_model_deltaslope

		*******************
		***				***
		*** RHO = -0.25  ***
		***				***
		*******************
	
		model -0.25				
		
		return scalar stay_model2_cons  	 	 = $stay_model_cons
		return scalar stay_model2_slope 	 	 = $stay_model_slope
		return scalar qual_model2_cons  	 	 = $qual_model_cons
		return scalar qual_model2_slope 	 	 = $qual_model_slope
		return scalar qual_model2_after 	 	 = $qual_model_after
		return scalar qual_model2_deltaslope  	 = $qual_model_deltaslope
		
		*******************
		***				***
		*** RHO = -0.1  ***
		***				***
		*******************
	
		model -0.1				
		
		return scalar stay_model3_cons  	 	 = $stay_model_cons
		return scalar stay_model3_slope 	 	 = $stay_model_slope
		return scalar qual_model3_cons  	 	 = $qual_model_cons
		return scalar qual_model3_slope 	 	 = $qual_model_slope
		return scalar qual_model3_after 	 	 = $qual_model_after
		return scalar qual_model3_deltaslope  	 = $qual_model_deltaslope
		
		*******************
		***				***
		*** RHO =  0    ***
		***				***
		*******************
	
		model 0				
		
		return scalar stay_model4_cons  	 	 = $stay_model_cons
		return scalar stay_model4_slope 	 	 = $stay_model_slope
		return scalar qual_model4_cons  	 	 = $qual_model_cons
		return scalar qual_model4_slope 	 	 = $qual_model_slope
		return scalar qual_model4_after 	 	 = $qual_model_after
		return scalar qual_model4_deltaslope  	 = $qual_model_deltaslope

		*******************
		***				***
		*** RHO =  0.1  ***
		***				***
		*******************
	
		model 0.1				
		
		return scalar stay_model5_cons  	 	 = $stay_model_cons
		return scalar stay_model5_slope 	 	 = $stay_model_slope
		return scalar qual_model5_cons  	 	 = $qual_model_cons
		return scalar qual_model5_slope 	 	 = $qual_model_slope
		return scalar qual_model5_after 	 	 = $qual_model_after
		return scalar qual_model5_deltaslope  	 = $qual_model_deltaslope
		
		*******************
		***				***
		*** RHO = 0.25  ***
		***				***
		*******************
	
		model 0.25				
		
		return scalar stay_model6_cons  	 	 = $stay_model_cons
		return scalar stay_model6_slope 	 	 = $stay_model_slope
		return scalar qual_model6_cons  	 	 = $qual_model_cons
		return scalar qual_model6_slope 	 	 = $qual_model_slope
		return scalar qual_model6_after 	 	 = $qual_model_after
		return scalar qual_model6_deltaslope  = $qual_model_deltaslope

		*******************
		***				***
		*** RHO = 0.5   ***
		***				***
		*******************
	
		model 0.5				
		
		return scalar stay_model7_cons  	 	 = $stay_model_cons
		return scalar stay_model7_slope 	 	 = $stay_model_slope
		return scalar qual_model7_cons  	 	 = $qual_model_cons
		return scalar qual_model7_slope 	 	 = $qual_model_slope
		return scalar qual_model7_after 	 	 = $qual_model_after
		return scalar qual_model7_deltaslope  = $qual_model_deltaslope

	end
	
		
	bootstrap 	stay_actual_cons =r(stay_actual_cons) ///
				stay_model1_cons=r(stay_model1_cons) ///	
				stay_model2_cons=r(stay_model2_cons) ///	
				stay_model3_cons=r(stay_model3_cons) ///	
				stay_model4_cons=r(stay_model4_cons) ///	
				stay_model5_cons=r(stay_model5_cons) ///	
				stay_model6_cons=r(stay_model6_cons) ///	
				stay_model7_cons=r(stay_model7_cons) ///	
				stay_actual_slope=r(stay_actual_slope) ///	
				stay_model1_slope=r(stay_model1_slope) ///	
				stay_model2_slope=r(stay_model2_slope) ///	
				stay_model3_slope=r(stay_model3_slope) ///	
				stay_model4_slope=r(stay_model4_slope) ///	
				stay_model5_slope=r(stay_model5_slope) ///	
				stay_model6_slope=r(stay_model6_slope) ///	
				stay_model7_slope=r(stay_model7_slope) ///	
				qual_actual_cons =r(qual_actual_cons) ///	
				qual_model1_cons=r(qual_model1_cons) ///	
				qual_model2_cons=r(qual_model2_cons) ///	
				qual_model3_cons=r(qual_model3_cons) ///	
				qual_model4_cons=r(qual_model4_cons) ///	
				qual_model5_cons=r(qual_model5_cons) ///	
				qual_model6_cons=r(qual_model6_cons) ///	
				qual_model7_cons=r(qual_model7_cons) ///	
				qual_actual_slope=r(qual_actual_slope) ///	
				qual_model1_slope=r(qual_model1_slope) ///	
				qual_model2_slope=r(qual_model2_slope) ///	
				qual_model3_slope=r(qual_model3_slope) ///	
				qual_model4_slope=r(qual_model4_slope) ///	
				qual_model5_slope=r(qual_model5_slope) ///	
				qual_model6_slope=r(qual_model6_slope) ///	
				qual_model7_slope=r(qual_model7_slope) ///	
				qual_actual_after=r(qual_actual_after) ///	
				qual_model1_after=r(qual_model1_after) ///	
				qual_model2_after=r(qual_model2_after) ///	
				qual_model3_after=r(qual_model3_after) ///	
				qual_model4_after=r(qual_model4_after) ///	
				qual_model5_after=r(qual_model5_after) ///	
				qual_model6_after=r(qual_model6_after) ///	
				qual_model7_after=r(qual_model7_after) ///	
				qual_actual_deltaslope=r(qual_actual_deltaslope) ///	
				qual_model1_deltaslope=r(qual_model1_deltaslope) ///		
				qual_model2_deltaslope=r(qual_model2_deltaslope) ///		
				qual_model3_deltaslope=r(qual_model3_deltaslope) ///		
				qual_model4_deltaslope=r(qual_model4_deltaslope) ///		
				qual_model5_deltaslope=r(qual_model5_deltaslope) ///		
				qual_model6_deltaslope=r(qual_model6_deltaslope) ///		
				qual_model7_deltaslope=r(qual_model7_deltaslope) ///		
				, reps(1000) seed(15356410) nodrop: boot_model

				forvalues i = 1/7 {
					
					foreach construct in cons slope after deltaslope {
						test (qual_model`i'_`construct' = qual_actual_`construct')
						local qual_`construct'`i' = r(p)
					}
					
					foreach construct in cons slope {
						test (stay_model`i'_`construct' = stay_actual_`construct')
						local stay_`construct'`i' = r(p)
					}
					
					test (qual_model`i'_cons  = qual_model`i'_after)
					local ROSLA_intercept`i' = r(p)
					
					test (qual_model`i'_slope = qual_model`i'_deltaslope)
					local ROSLA_slope`i' = r(p)
					
				}
